subroutine Possion3DH201
    use Mesh3DH20_m
    use SMatrix_m
    use FES3DH201_m
    use IterSSolver_m
    implicit none
    type(mesh3DH8) :: mesh3D0
    type(mesh3DH20) :: mesh3D
    type(FES3DH201) :: fs
    type(SMatrix) :: A
    type(IterSSolver) :: solver

    real*8, allocatable :: x(:), RHS(:)

    integer :: label(6), i

    integer :: iter

    call mesh3D0.Cube(20, 20, 20)
    call mesh3D.initMesh(mesh3D0)
    call fs.initFES(mesh3D)

!    call print3DH20DOFS(fs)
!    call print3DH20Elememts(fs)
!    call print3DH20Boundarys(fs)

    call A.initSMatrix(fs%nDof)
    write(*, *) "The Number of DOF", fs%nDof
    allocate(RHS(fs%nDof), x(fs%nDof))
    
    RHS = 0.0D0
    x = 0.0D0
    label = (/1, 2, 3, 4, 5, 6/)
    call fs.assembleHeatStiff(constFunc3D, 0.0D0, COMMON_2D, A)
!    call sortSparseMatrix(A)
!    call displaySparseMatrix(A)

    call fs.assembleHeatSource(source3D, 0.0D0, COMMON_2D, RHS)

!    write(*, *)"RHS"
!    write(*, "(F20.12)")RHS

    write(*, *)"Before Apply BC"
    write(*, *) "The Size Of SparseMatrix", A%NNZ

   call fs.applyBCZero(A, RHS, label)

    write(*, *)"After Apply BC"
    write(*, *) "The Size Of SparseMatrix", A%NNZ

    call solver.initIterSSolver
    call solver.PCGSSOR(A, RHS, x, 1.5D0, 2)

    write(*, *) "Solution:"
    write(*, *) "min : ", minval(x), "max: ", maxval(x)
    
    !call writeToTxt("Q20ForMatlab.txt", len("Q20ForMatlab.txt"), fs, x)
!   call exportMeshSolution("solu3DH.dat", len("solu3DH.dat"), mesh3D, x)
!    write(*, *) GaussPoints3
    call deallocateSMatrix(A)
    deallocate(RHS, x)
contains

real*8 function constFunc3D(xyz, nRegion, time)
    implicit none
    real*8, intent(in) :: xyz(3), time
    integer, intent(in) :: nRegion
    constFunc3D = 1.0D0
end function

real*8 function source3D(xyz, nRegion, time)
    implicit none
    real*8, intent(in) :: xyz(3), time
    integer, intent(in) :: nRegion
    
    real*8 :: x, y, z
    x = xyz(1); y = xyz(2); z = xyz(3)
    source3D = 2*x*(1-x)*y*(1-y) + 2*y*(1-y)*z*(1-z) + 2*x*(1-x)*z*(1-z)
end function

end subroutine